# R source code for the course project of STATS 240P
# Set up the experiment directory
rm(list=ls())

library(ggplot2)
library(xts)
library(fGarch)
source("Forecast.R")
source("BL.R") # Black-Litterman
source("MV.R") # Markowitz Optimization
source("Eval.R") # Evaluation

FF6 <- read.table("../data/FF6Portfolios.txt", header=T, skip=3)
are.log <- TRUE
if (are.log) {
	FF6$smlo_vwret <- log(FF6$smlo_vwret + 1)
	FF6$smme_vwret <- log(FF6$smme_vwret + 1)
	FF6$smhi_vwret <- log(FF6$smhi_vwret + 1)
	FF6$bilo_vwret <- log(FF6$bilo_vwret + 1)
	FF6$bime_vwret <- log(FF6$bime_vwret + 1)
	FF6$bihi_vwret <- log(FF6$bihi_vwret + 1)
}

NN <- nrow(FF6) # Set this as a variable for quick experiments.
#NN <- 66
use.mvp <- TRUE
ff.test <- FF6[61:NN,]
mv.plugin.intro() #Charts in section 1; uncomment if need to re-run.

# Experiment 1: using expanding window + Markowitz Plug-in Estimates.
mv.plugin <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 1: Markowitz Plug-in Estimates using Expanding Window")
for (t0 in 61:NN) {
	ff.train <- FF6[1:(t0-1),] # use the expanding window for training
	ret <- data.frame(ff.train$smlo_vwret, ff.train$smme_vwret, ff.train$smhi_vwret, ff.train$bilo_vwret, ff.train$bime_vwret, ff.train$bihi_vwret)
	names(ret) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")
	mu <- as.matrix(mean(ret))
	sigma <- cov(ret)
	aversion <- compute.risk.aversion()
	
	w <- optimal.weights(mu, sigma, aversion)
	if (use.mvp) {
		w <- weights.mvp(mu, sigma)
	}
	mv.plugin[t0-60, 2:7] <- c(w)
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.transition.map(mv.plugin[,2:7], "../results/mv_transition.png")
	plot.abs.transition.map(mv.plugin[,2:7], "../results/mv_abs_transition.png")
} else {
	plot.transition.map(mv.plugin[,2:7], "../results/mv_transition_aversion.png")
	plot.abs.transition.map(mv.plugin[,2:7], "../results/mv_abs_transition_aversion.png")
}
return.mv.expand <- compute.returns(mv.plugin, ff.test, are.log)
sharpe.mv.expand <- get.sharpe.ratio(return.mv.expand)

# End of Experiment 1.

###  Black-Litterman Model ###

# Experiment 2: ARMA model + Expanding window
# We use an expanding window as training data, starting from 5 years of data.
# alloc.arima.expand: track weight allocation across test period
alloc.arima.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.arima.expand: track return prediction across test period
pred.arima.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.arima.expand: track standard error across test period
se.arima.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 2: B-L + ARMA model using Expanding Window")
for (t0 in 61:NN) {
	ff.train <- FF6[1:(t0-1),] # use the expanding window for training
	# ARMA
	pred <- forecast.arima(ff.train)
	Q <- pred[["Q"]]
	pred.arima.expand[t0-60, 2:7] <- Q
	SE <- pred[["SE"]]
	se.arima.expand[t0-60, 2:7] <- SE
	
	# Black-litterman
	P <- diag(6)
	bl.arima <- black.litterman(ff.train, P, as.matrix(Q)) 
	mu.bl.arima <- bl.arima[["mu"]]
	sigma.bl.arima <- bl.arima[["sigma"]]
	aversion <- bl.arima[["aversion"]]
	alloc.arima.expand[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.arima, sigma.bl.arima, aversion))
	if (use.mvp) {
		alloc.arima.expand[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.arima, sigma.bl.arima))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.arima.expand, pred.arima.expand, se.arima.expand, ff.test, "arma_expand")
} else {
	plot.performance(alloc.arima.expand, pred.arima.expand, se.arima.expand, ff.test, "arma_expand_aversion")
}
return.arma.expand <- compute.returns(alloc.arima.expand, ff.test, are.log)
sharpe.arma.expand <- get.sharpe.ratio(return.arma.expand)
# End of Experiment 2

# Experiment 3: AR(1) + Garch(1,1) model + Expanding window
# We use an expanding window as training data, starting from 5 years of data.
# alloc.garch.expand: track weight allocation across test period
alloc.garch.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.garch.expand: track return prediction across test period
pred.garch.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.garch.expand: track standard error across test period
se.garch.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 3: B-L + AR(1) + Garch(1,1) model using Expanding Window")
for (t0 in 61:NN) {
	ff.train <- FF6[1:(t0-1),] # use the expanding window for training	
	ret <- data.frame(ff.train$smlo_vwret, ff.train$smme_vwret, ff.train$smhi_vwret, ff.train$bilo_vwret, ff.train$bime_vwret, ff.train$bihi_vwret)
	names(ret) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")

	# Build AR(1)-GARCH(1,1) models
	smlo.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),1], trace=FALSE)
	smme.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),2], trace=FALSE)
	smhi.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),3], trace=FALSE)
	bilo.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),4], trace=FALSE)
	bime.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),5], trace=FALSE)
	bihi.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[1:(t0-1),6], trace=FALSE)

	# Generate AR-GARCH prediction and load into view matrix, Q (6 x 1 column vector)
	Q<-matrix(0,6,1)
	SE <- rep(NA, 6)
	Q[1,1]<-predict(smlo.model,n.ahead=1)$meanForecast
	SE[1]<-predict(smlo.model,n.ahead=1)$standardDeviation
	Q[2,1]<-predict(smme.model,n.ahead=1)$meanForecast
	SE[2]<-predict(smme.model,n.ahead=1)$standardDeviation
	Q[3,1]<-predict(smhi.model,n.ahead=1)$meanForecast
	SE[3]<-predict(smhi.model,n.ahead=1)$standardDeviation
	Q[4,1]<-predict(bilo.model,n.ahead=1)$meanForecast
	SE[4]<-predict(bilo.model,n.ahead=1)$standardDeviation
	Q[5,1]<-predict(bime.model,n.ahead=1)$meanForecast
	SE[5]<-predict(bime.model,n.ahead=1)$standardDeviation
	Q[6,1]<-predict(bihi.model,n.ahead=1)$meanForecast
	SE[6]<-predict(bihi.model,n.ahead=1)$standardDeviation
	
	pred.garch.expand[t0-60, 2:7] <- as.vector(Q)
	se.garch.expand[t0-60, 2:7] <- SE
	
	# Black-litterman
	P <- diag(6)
	bl.garch <- black.litterman(ff.train, P, Q) 
	mu.bl.garch <- bl.garch[["mu"]]
	sigma.bl.garch <- bl.garch[["sigma"]]
	aversion <- bl.garch[["aversion"]]
	
	alloc.garch.expand[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.garch, sigma.bl.garch, aversion))
	if (use.mvp) {
		alloc.garch.expand[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.garch, sigma.bl.garch))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.garch.expand, pred.garch.expand, se.garch.expand, ff.test, "garch_expand")
} else {
	plot.performance(alloc.garch.expand, pred.garch.expand, se.garch.expand, ff.test, "garch_expand_aversion")
}
return.garch.expand <- compute.returns(alloc.garch.expand, ff.test, are.log)
sharpe.garch.expand <- get.sharpe.ratio(return.garch.expand)
# End of Experiment 3

# Experiment 4: VIX model + Expanding window
# We use an expanding window as training data, starting from 5 years of data.
#
# 1) Load change in VIX data (dVIX)
#
# 2) Convert log-returns of FF portfolios and change in VIX data to time series data
dVIX <- read.csv("../data/vix.csv", header=T)
FFRet <- data.frame(FF6$date,FF6$smlo_vwret,FF6$smme_vwret,FF6$smhi_vwret, FF6$bilo_vwret, FF6$bime_vwret, FF6$bihi_vwret)
FFRet.ts<-ts(FFRet[,2:7],start=c(1992,1),freq=12) # assuming FFRet[,2:7] holds log returns for all 6 FF portfolios
dVIX.ts <-ts(dVIX,start=c(1992,1),freq=12)
# alloc.vix.expand: track weight allocation across test period
alloc.vix.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.vix.expand: track return prediction across test period
pred.vix.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.vix.expand: track standard error across test period
se.vix.expand <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 4: B-L + VIX model using Expanding Window")
for (t0 in 61:NN) {
 	# 3) Build single factor linear model that describes FF portfolio movement by the change in VIX
 	vix.train <- data.frame(x=dVIX.ts[1:(t0-1)], y1=FFRet.ts[1:(t0-1),1], y2=FFRet.ts[1:(t0-1),2], y3=FFRet.ts[1:(t0-1),3], y4=FFRet.ts[1:(t0-1),4], y5=FFRet.ts[1:(t0-1),5], y6=FFRet.ts[1:(t0-1),6])
	SMLO.VIX <- lm(y1 ~ x, data=vix.train)
	SMME.VIX <- lm(y2 ~ x, data=vix.train)
	SMHI.VIX <- lm(y3 ~ x, data=vix.train)
	BILO.VIX <- lm(y4 ~ x, data=vix.train)
	BIME.VIX <- lm(y5 ~ x, data=vix.train)
	BIHI.VIX <- lm(y6 ~ x, data=vix.train)
	
	# VIX model
	
	# 4) Build ARIMA model for change in VIX
	
	# Identify and select best ARIMA model for VIX using Cai's code to find minimum AIC
	pq <- arma.model.selection(dVIX.ts[1:(t0-1)], 3, 3, 0)
	VIX.ARIMA <-arima(dVIX.ts[1:(t0-1)],order=c(pq[["p"]],0,pq[["q"]]))

	# 5) Generate t+1 VIX prediction to predict FF portfolio next time step return and load into view matrix, Q (6 x 1 column vector)
	Q<-matrix(0,6,1)
	SE <- rep(NA, 6)
	new<-data.frame(x<-predict(VIX.ARIMA,n.ahead=1)$pred) # t+1 VIX prediction
	Q[1,1]<-predict(SMLO.VIX,newdata=new, se.fit=TRUE)$fit
	SE[1]<-predict(SMLO.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[2,1]<-predict(SMME.VIX,newdata=new, se.fit=TRUE)$fit
	SE[2]<-predict(SMME.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[3,1]<-predict(SMHI.VIX,newdata=new, se.fit=TRUE)$fit
	SE[3]<-predict(SMHI.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[4,1]<-predict(BILO.VIX,newdata=new, se.fit=TRUE)$fit
	SE[4]<-predict(BILO.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[5,1]<-predict(BIME.VIX,newdata=new, se.fit=TRUE)$fit
	SE[5]<-predict(BIME.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[6,1]<-predict(BIHI.VIX,newdata=new, se.fit=TRUE)$fit
	SE[6]<-predict(BIHI.VIX,newdata=new, se.fit=TRUE)$se.fit
	pred.vix.expand[t0-60, 2:7] <- as.vector(Q)
	se.vix.expand[t0-60, 2:7] <- SE
	# Black-litterman
	P <- diag(6)
	bl.vix <- black.litterman(ff.train, P, Q) 
	mu.bl.vix <- bl.vix[["mu"]]
	sigma.bl.vix <- bl.vix[["sigma"]]
	aversion <- bl.vix[["aversion"]]
	alloc.vix.expand[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.vix, sigma.bl.vix, aversion))
	if (use.mvp) {
		alloc.vix.expand[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.vix, sigma.bl.vix))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.vix.expand, pred.vix.expand, se.vix.expand, ff.test, "vix_expand")
} else {
	plot.performance(alloc.vix.expand, pred.vix.expand, se.vix.expand, ff.test, "vix_expand_aversion")
}
return.vix.expand <- compute.returns(alloc.vix.expand, ff.test, are.log=FALSE)
sharpe.vix.expand <- get.sharpe.ratio(return.vix.expand)

# End of Experiment 4

# Experiment 5: ARMA model + Rolling window
# We use an rolling window as training data, starting from 5 years of data.
# alloc.arima.roll: track weight allocation across test period
alloc.arima.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.arima.roll: track return prediction across test period
pred.arima.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.arima.roll: track standard error across test period
se.arima.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 5: B-L + ARMA model using Rolling Window")
for (t0 in 61:NN) {
	ff.train <- FF6[(t0-60):(t0-1),] # use the expanding window for training
	# ARMA
	pred <- forecast.arima(ff.train)
	Q <- pred[["Q"]]
	pred.arima.roll[t0-60, 2:7] <- Q
	SE <- pred[["SE"]]
	se.arima.roll[t0-60, 2:7] <- SE
	
	# Black-litterman
	P <- diag(6)
	bl.arima <- black.litterman(ff.train, P, as.matrix(Q)) 
	mu.bl.arima <- bl.arima[["mu"]]
	sigma.bl.arima <- bl.arima[["sigma"]]
	aversion <- bl.arima[["aversion"]]
	alloc.arima.roll[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.arima, sigma.bl.arima, aversion))
	if (use.mvp) {
		alloc.arima.roll[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.arima, sigma.bl.arima))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.arima.roll, pred.arima.roll, se.arima.roll, ff.test, "arma_roll")
} else {
	plot.performance(alloc.arima.roll, pred.arima.roll, se.arima.roll, ff.test, "arma_roll_aversion")
}
return.arma.roll <- compute.returns(alloc.arima.roll, ff.test, are.log=FALSE)
sharpe.arma.roll <- get.sharpe.ratio(return.arma.roll)
# End of Experiment 5

# Experiment 6: AR(1) + Garch(1,1) model + Rolling window
# We use an rolling window as training data, starting from 5 years of data.
# alloc.garch.roll: track weight allocation across test period
alloc.garch.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.garch.roll: track return prediction across test period
pred.garch.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.garch.roll: track standard error across test period
se.garch.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 6: B-L + AR(1) + Garch(1,1) model using Rolling Window")
for (t0 in 61:NN) {
	ff.train <- FF6[(t0-60):(t0-1),] # use the expanding window for training	
	ret <- data.frame(ff.train$smlo_vwret, ff.train$smme_vwret, ff.train$smhi_vwret, ff.train$bilo_vwret, ff.train$bime_vwret, ff.train$bihi_vwret)
	names(ret) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")

	# Build AR(1)-GARCH(1,1) models
	smlo.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,1], trace=FALSE)
	smme.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,2], trace=FALSE)
	smhi.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,3], trace=FALSE)
	bilo.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,4], trace=FALSE)
	bime.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,5], trace=FALSE)
	bihi.model <- garchFit(formula=as.formula("~arma(1,0)+garch(1,1)"),data=ret[,6], trace=FALSE)

	# Generate AR-GARCH prediction and load into view matrix, Q (6 x 1 column vector)
	Q<-matrix(0,6,1)
	SE <- rep(NA, 6)
	Q[1,1]<-predict(smlo.model,n.ahead=1)$meanForecast
	SE[1]<-predict(smlo.model,n.ahead=1)$standardDeviation
	Q[2,1]<-predict(smme.model,n.ahead=1)$meanForecast
	SE[2]<-predict(smme.model,n.ahead=1)$standardDeviation
	Q[3,1]<-predict(smhi.model,n.ahead=1)$meanForecast
	SE[3]<-predict(smhi.model,n.ahead=1)$standardDeviation
	Q[4,1]<-predict(bilo.model,n.ahead=1)$meanForecast
	SE[4]<-predict(bilo.model,n.ahead=1)$standardDeviation
	Q[5,1]<-predict(bime.model,n.ahead=1)$meanForecast
	SE[5]<-predict(bime.model,n.ahead=1)$standardDeviation
	Q[6,1]<-predict(bihi.model,n.ahead=1)$meanForecast
	SE[6]<-predict(bihi.model,n.ahead=1)$standardDeviation
	
	pred.garch.roll[t0-60, 2:7] <- as.vector(Q)
	se.garch.roll[t0-60, 2:7] <- SE
	
	# Black-litterman
	P <- diag(6)
	bl.garch <- black.litterman(ff.train, P, Q) 
	mu.bl.garch <- bl.garch[["mu"]]
	sigma.bl.garch <- bl.garch[["sigma"]]
	aversion <- bl.garch[["aversion"]]
	alloc.garch.roll[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.garch, sigma.bl.garch, aversion))
	if (use.mvp) {
		alloc.garch.roll[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.garch, sigma.bl.garch))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.garch.roll, pred.garch.roll, se.garch.roll, ff.test, "garch_roll")
} else {
	plot.performance(alloc.garch.roll, pred.garch.roll, se.garch.roll, ff.test, "garch_roll_aversion")
}
return.garch.roll <- compute.returns(alloc.garch.roll, ff.test, are.log=FALSE)
sharpe.garch.roll <- get.sharpe.ratio(return.garch.roll)
# End of Experiment 6

# Experiment 7: VIX model + Rolling window
# We use an rolling window as training data, starting from 5 years of data.
#
# 1) Load change in VIX data (dVIX)
#
# 2) Convert log-returns of FF portfolios and change in VIX data to time series data
dVIX <- read.csv("../data/vix.csv", header=T)
FFRet <- data.frame(FF6$date,FF6$smlo_vwret,FF6$smme_vwret,FF6$smhi_vwret, FF6$bilo_vwret, FF6$bime_vwret, FF6$bihi_vwret)
FFRet.ts<-ts(FFRet[,2:7],start=c(1992,1),freq=12) # assuming FFRet[,2:7] holds log returns for all 6 FF portfolios
dVIX.ts <-ts(dVIX,start=c(1992,1),freq=12)
# alloc.vix.roll: track weight allocation across test period
alloc.vix.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# pred.vix.roll: track return prediction across test period
pred.vix.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
# se.vix.roll: track standard error across test period
se.vix.roll <- data.frame(date=FF6$date[61:NN], smlo=rep(NA, NN-60), smme=rep(NA, NN-60), smhi=rep(NA, NN-60), bilo=rep(NA, NN-60), bime=rep(NA, NN-60), bihi=rep(NA, NN-60))
print("Running experiment 7: B-L + VIX model using Rolling Window")
for (t0 in 61:NN) {
 	# 3) Build single factor linear model that describes FF portfolio movement by the change in VIX
 	vix.train <- data.frame(x=dVIX.ts[(t0-60):(t0-1)], y1=FFRet.ts[(t0-60):(t0-1),1], y2=FFRet.ts[(t0-60):(t0-1),2], y3=FFRet.ts[(t0-60):(t0-1),3], y4=FFRet.ts[(t0-60):(t0-1),4], y5=FFRet.ts[(t0-60):(t0-1),5], 	  y6=FFRet.ts[(t0-60):(t0-1),6])
	SMLO.VIX <- lm(y1 ~ x, data=vix.train)
	SMME.VIX <- lm(y2 ~ x, data=vix.train)
	SMHI.VIX <- lm(y3 ~ x, data=vix.train)
	BILO.VIX <- lm(y4 ~ x, data=vix.train)
	BIME.VIX <- lm(y5 ~ x, data=vix.train)
	BIHI.VIX <- lm(y6 ~ x, data=vix.train)
	
	# VIX model
	
	# 4) Build ARIMA model for change in VIX
	
	# Identify and select best ARIMA model for VIX using Cai's code to find minimum AIC
	pq <- arma.model.selection(dVIX.ts[(t0-60):(t0-1)], 3, 3, 0)
	VIX.ARIMA <-arima(dVIX.ts[(t0-60):(t0-1)],order=c(pq[["p"]],0,pq[["q"]]), method="ML")

	# 5) Generate t+1 VIX prediction to predict FF portfolio next time step return and load into view matrix, Q (6 x 1 column vector)
	Q<-matrix(0,6,1)
	SE <- rep(NA, 6)
	new<-data.frame(x<-predict(VIX.ARIMA,n.ahead=1)$pred) # t+1 VIX prediction
	Q[1,1]<-predict(SMLO.VIX,newdata=new, se.fit=TRUE)$fit
	SE[1]<-predict(SMLO.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[2,1]<-predict(SMME.VIX,newdata=new, se.fit=TRUE)$fit
	SE[2]<-predict(SMME.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[3,1]<-predict(SMHI.VIX,newdata=new, se.fit=TRUE)$fit
	SE[3]<-predict(SMHI.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[4,1]<-predict(BILO.VIX,newdata=new, se.fit=TRUE)$fit
	SE[4]<-predict(BILO.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[5,1]<-predict(BIME.VIX,newdata=new, se.fit=TRUE)$fit
	SE[5]<-predict(BIME.VIX,newdata=new, se.fit=TRUE)$se.fit
	Q[6,1]<-predict(BIHI.VIX,newdata=new, se.fit=TRUE)$fit
	SE[6]<-predict(BIHI.VIX,newdata=new, se.fit=TRUE)$se.fit
	pred.vix.roll[t0-60, 2:7] <- as.vector(Q)
	se.vix.roll[t0-60, 2:7] <- SE
	# Black-litterman
	P <- diag(6)
	bl.vix <- black.litterman(ff.train, P, Q) 
	mu.bl.vix <- bl.vix[["mu"]]
	sigma.bl.vix <- bl.vix[["sigma"]]
	aversion <- bl.vix[["aversion"]]
	alloc.vix.roll[t0-60, 2:7] <- unlist(optimal.weights(mu.bl.vix, sigma.bl.vix, aversion))
	if (use.mvp) {
		alloc.vix.roll[t0-60, 2:7] <- unlist(weights.mvp(mu.bl.vix, sigma.bl.vix))
	}
	print(FF6[t0, 1])
}

if (use.mvp) {
	plot.performance(alloc.vix.roll, pred.vix.roll, se.vix.roll, ff.test, "vix_roll")
} else {
	plot.performance(alloc.vix.roll, pred.vix.roll, se.vix.roll, ff.test, "vix_roll_aversion")
}
return.vix.roll <- compute.returns(alloc.vix.roll, ff.test, are.log=FALSE)
sharpe.vix.roll <- get.sharpe.ratio(return.vix.roll)

# End of Experiment 7

# Final plots for results using two metrics:
# Cumulative returns:
#   return.mvp, return.arma
returns <- data.frame(date=ff.test$date, mv_expand=return.mv.expand, arma_expand=return.arma.expand, garch_expand=return.garch.expand, vix_expand=return.vix.expand,
	arma_roll=return.arma.roll, garch_roll=return.garch.roll, vix_roll=return.vix.roll)
if (use.mvp) {
	png("../results/eval_return.png")
} else {
	png("../results/eval_return_aversion.png")
}
ggplot(returns, aes(x=as.Date(as.character(date), format="%Y%m%d"))) + 
	geom_line(aes(y=mv_expand, colour="mv-plugin-expand")) + 
	geom_line(aes(y=arma_expand, colour="bl-arma-expand")) +
	geom_line(aes(y=garch_expand, colour="bl-garch-expand")) + 
	geom_line(aes(y=vix_expand, colour="bl-vix-expand")) +
	geom_line(aes(y=arma_roll, colour="bl-arma-roll")) +
	geom_line(aes(y=garch_roll, colour="bl-garch-roll")) + 
	geom_line(aes(y=vix_roll, colour="bl-vix-roll")) +
	scale_x_date() + xlab("date") + ylab("cumulative return")

dev.off()

sharpes <- data.frame(date=ff.test$date, mv_expand=sharpe.mv.expand, arma_expand=sharpe.arma.expand, garch_expand=sharpe.garch.expand, vix_expand=sharpe.vix.expand,
arma_roll=sharpe.arma.roll, garch_roll=sharpe.garch.roll, vix_roll=sharpe.vix.roll)
if (use.mvp) {
	png("../results/eval_sharpe.png")
} else {
	png("../results/eval_sharpe_aversion.png")
}
ggplot(sharpes, aes(x=as.Date(as.character(date), format="%Y%m%d"))) + 
	geom_line(aes(y=mv_expand, colour="mv-plugin-expand")) + 
	geom_line(aes(y=arma_expand, colour="bl-arma-expand")) + 
	geom_line(aes(y=garch_expand, colour="bl-garch-expand")) +
	geom_line(aes(y=vix_expand, colour="bl-vix-expand")) +
	geom_line(aes(y=arma_roll, colour="bl-arma-roll")) + 
	geom_line(aes(y=garch_roll, colour="bl-garch-roll")) +
	geom_line(aes(y=vix_roll, colour="bl-vix-roll")) +
	scale_x_date() + xlab("date") + ylab("sharpe ratio")
dev.off()
